clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Luca, how many years are used for the initilalization? We should use 10
%%% years, 1948-1958 as in CPS

%***********PARAMETERS SPECIFICATION Time Varying Var**************
gbit = 1;
maxgbitInit     =  12000;        % total number of draws
burnInit        =  10000;        % burn in period

maxgbit         =   2500;        % total number of draws
burn            =    500;        % burn in period

jump            =     5;           % collect one each jump draws


horizon = 12;
p = 2;


tight_Q = 0.00001;
tight_W = 0.00001;
tight_S = 0.001;

post_index = 1;
trs=0;
trsMax=10;
MaxEig=1.0;
MaxRoot=MaxEig+1;
c=1;
L=p;



%%******************************************************************
%%% Loading the realtime database for the unemployment rate (quarterly vintages monthly frequency)
[A B] = xlsread('data_us.xls');

INFL =(log(A(2:end,2))-log(A(1:end-1,2)))*400;
IR =(log(A(2:end,3))-log(A(1:end-1,3)))*400;
M2 =(log(A(2:end,4))-log(A(1:end-1,4)))*400;
GDP =(log(A(2:end,5))-log(A(1:end-1,5)))*400;


qq =datevec(A(1:end,1));
dd=find(qq(:,1)==0);
date = [datevec(datenum('30-Dec-1799')+A(2:dd-1,1)); datevec(datenum('30-Dec-1899')+A(dd:end,1))];
DATE = [date(:,1) date(:,2)];
TT = size(DATE,1);

DATA = INFL;


%******************************************************************


start_eval_y  = 1900;   start_eval_m = 3;  %% start he evaluation from 1970Q1
start_sim = find(DATE(:,1) == start_eval_y & DATE(:,2) == start_eval_m)-1;

T =size(DATA,1);
EstWind = 80;
Z = [INFL M2 GDP];



EstWind = 50;

tic
count = 0;
for j = 152:T-horizon
    j
    count = count +1;
    disp('--------------------------')
    disp('--------------------------')
    disp(['Now running the forecast at ',num2str(DATE(j,1)),'Q',num2str(DATE(j,2))])



    infl = INFL(1:j,:);
    m2 = M2(1:j,:);
    gdp = GDP(1:j,:);

    X = [infl m2 gdp];

    PRED_var_rec = VAR_fore(X,p,horizon);
    PRED_var_rol = VAR_fore(X(end-EstWind+1:end,:),p,horizon);


    PRED_ar_rec(1,:) = VAR_fore(X(:,1),p,horizon);
    PRED_ar_rol(1,:) = VAR_fore(X(end-EstWind+1:end,1),p,horizon);




    MixingKSC
    %clear PRED_tv_var PB
    if count == 1
        %% girare con init, chiedere parametri come output, var resid very
        %% important since it is going to be the input for the next start
        [PRED_tv_var, PB,prevR, prevy2star,prevB,prevA] = ...
            TvVarEstKSC(X,tight_Q,tight_W,tight_S,gbit,maxgbitInit,burnInit,jump,post_index,trs,trsMax,MaxEig,MaxRoot,c,L,horizon,ksc);

    else

        [PRED_tv_var, PB, prevR, prevy2star,prevB,prevA] = ...
            TvVarInputEstKSC(X,tight_Q,tight_W,tight_S,gbit,maxgbit,burn,jump,post_index,trs,trsMax,MaxEig,MaxRoot,c,L,horizon,prevR,prevy2star,prevB,prevA,ksc);
   
    end;

    for jg=1:size(PB,3)
        MatA=inv(eye(L*size(X,2))-companion(PB(:,end,jg),L,size(X,2),c));
        LM(:,jg)=MatA(1:size(X,2),1:size(X,2))*PB(1:L*size(X,2)+c:end,end,jg);
    end

    LM_d(j,:) = prctile(squeeze(LM)',16);
    LM_m(j,:) = prctile(squeeze(LM)',50);
    LM_u(j,:) = prctile(squeeze(LM)',84);




    clear PRED_tv_ar


    if count == 1
        [PRED_tv_ar, PBUni,prevRUni, prevy2starUni,prevBUni,prevAUni] = ...
            TvVarEstKSC(X(:,1),tight_Q,tight_W,tight_S,gbit,maxgbitInit,burnInit,jump,post_index,trs,trsMax,MaxEig,MaxRoot,c,L,horizon,ksc);
    else

           [PRED_tv_ar, PBUni, prevRUni, prevy2starUni,prevBUni,prevAUni] = ...
            TvVarInputEstKSC(X(:,1),tight_Q,tight_W,tight_S,gbit,maxgbit,burn,jump,post_index,trs,trsMax,MaxEig,MaxRoot,c,L,horizon,prevRUni,prevy2starUni,prevBUni,prevAUni,ksc);
           
    end;

    for jg=1:size(PBUni,3)
        LMUni(:,jg)=(1-ones(1,L)*PBUni(c+1:end,end,jg))^(-1)*PBUni(1,end,jg);
    end

    LM_uni_d(j,1) = prctile(squeeze(LMUni)',16);
    LM_uni_m(j,1) = prctile(squeeze(LMUni)',50);
    LM_uni_u(j,1) = prctile(squeeze(LMUni)',84);

    
    for h = 1:horizon

        fore_var_rec_cum{h}(j+h,:)    = mean(PRED_var_rec(:,1:h),2)';
        fore_var_rol_cum{h}(j+h,:)    = mean(PRED_var_rol(:,1:h),2)';

        fore_ar_rec_cum{h}(j+h,:)    = mean(PRED_ar_rec(:,1:h),2)';
        fore_ar_rol_cum{h}(j+h,:)    = mean(PRED_ar_rol(:,1:h),2)';

        fore_var_tv_d_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_var(:,1:h,:),2))',16);
        fore_var_tv_m_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_var(:,1:h,:),2))',50);
        fore_var_tv_u_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_var(:,1:h,:),2))',84);

        fore_ar_tv_d_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_ar(:,1:h,:),2))',16);
        fore_ar_tv_m_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_ar(:,1:h,:),2))',50);
        fore_ar_tv_u_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_ar(:,1:h,:),2))',84);

        fore_var_tv_m_cum_gibb{h} = squeeze(mean(PRED_tv_var(:,1:h,:),2));
        fore_ar_tv_m_cum_gibb{h} = squeeze(mean(PRED_tv_ar(:,1:h,:),2));
        
        fore_rw{h}(j+h,1)= mean(X(end-3:end,1));
        true_cum{h}(j+h,1) = mean(INFL(j+1:j+h,1));


    end

      eval(['save ', 'results_var_infl_m2_gdp_and_ar_',num2str(j), ' fore_var_tv_m_cum_gibb fore_ar_tv_m_cum_gibb fore_rw true_cum fore_var_rec_cum fore_var_rol_cum fore_ar_rec_cum fore_ar_rol_cum'])



end
%        eval(['save ', 'results_var_gdp_m2_and_ar_',num2str(j),...
%            'fore_var_tv_m_cum_gibb  fore_ar_tv_m_cum_gibb fore_rw{h} true_cum fore_var_rec_cum fore_var_rol_cum fore_ar_rec_cum fore_ar_rol_cum'])
%  save results_hist_infl_m2_prova1


